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Abstract. To support physically faithful simulation, numerical methods for hyperbolic conserva- 
tion laws are needed that efficiently mimic the constraints satisfied by exact solutions, including ma- 
terial conservation and positivity, while also maintaining high-order accuracy and numerical stability. 
Finite volume methods such as discontinuous Galerkin (DG) and weighted essentially non-oscillatory 
(WENO) schemes allow efficient high-order accuracy while maintaining conservation. Positivity lim- 
iters developed by Zhang and Shu and summarized in [Proc. R. Soc. A 467, 2752 (20f f)] ensure a 
minimum time step for which positivity of cell average quantities is maintained without sacrificing 
conservation or formal accuracy; this is achieved by linearly damping the deviation from the cell 
average just enough to enforce a cell positivity condition that requires positivity at boundary nodes 
and strategically chosen interior points. 

We assume that the set of positive states is convex; it follows that positivity is equivalent to scalar 
positivity of a collection of affine functionals. Based on this observation, we generalize the technique 
of Zhang and Shu to a framework that we call outflow positivity limiting: First, enforce positivity at 
boundary nodes. If wave speed desingularization is needed, cap wave speeds at physically justified 
maxima by using remapped states to calculate fluxes. Second, apply linear damping again to cap the 
boundary average of all positivity functionals at the maximum possible (relative to the cell average) 
for a scalar-valued representation positive in each mesh cell. This be done by enforcing positivity 
of the retentional, an afflne combination of the cell average and the boundary average, in the same 
way that Zhang and Shu would enforce positivity at a single point (and with similar computational 
expense). Third, limit the time step so that cell outflow is less than the initial cell content. 

We show that enforcing positivity at the interior points in Zhang and Shu's method is actually 
a means of capping boundary averages at the maximum possible for a positive solution. Capping 
boundary averages allows computational interior points to be chosen without sacrificing the guaran- 
teed positivity-preserving time step so as to optimize stabilization benefits relative to computational 
expense, e.g. by choosing points that coincide with nodal points of a DG scheme. 
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1. Introduction. Consider a variable-coefficient hyperbolic scalar partial differ- 
ential equation (PDE) of the form 

dtu{t,Ti) + V-f{t,x,u) = for xefi, i6M>o, (1.1) 

where c is the physical domain, t is the time coordinate, V • represents the 
divergence with respect to x, m : M>o x M.^ -> R is a conserved scalar quantity, and the 
flux function f : assumed differentiable. Assume that solutions 



of (1.1 1 are positivity-preserving: that is, if ■ti(0,x) > 0, then u(t,x) > for all t> 0. 
The context of this work will mostly be restricted to a single mesh cell in the interior 
of the domain, so we do not concern ourselves with boundary conditions for the PDE. 



A numerical method for equation (1.1) evolves a numerical solution U that ap- 



proximates u and is represented by a finite number of degrees of freedom; we say that 
the method is accurate if the error ||i7 - u||oo vanishes as the degrees of freedom are 
increased in an appropriate way. 
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Assumed properties and requirements 

We assume a hyperbolic conservation law of the form 

9tu(t,x) + V-f(t,x,u) = 0, (1.2) 

where hyperbolicity means that the flux Jacobian matrix (?, j) i-^ n- is diagonal- 
izable with real eigenvalues and a full set of eigenvectors for any unit direction vector 
n 6 E-°, where ||n|| = 1. Each eigenvalue is a wave speed. Valid physical solutions 



form of conservation law ( 1.2 1 



violate the differential form (1.2) at shock discontinuities but still satisfy the integral 



dt f u+ (f n-f = 0, (1.3) 

Jk JdK 

for arbitrary K. We assume that there is a convex set V (satisfying sV + {1- s)V c 
T-" Vs 6 [0, 1]) of positive states such that V is an invariant domain: if the initial data 
m(0, x) is in V then the physical solution remains in V for all time: u{t, x) e P Vt > 0. 
Furthermore, we assume that f ^ on the boundary of V. Positivity limiters yield 
a numerical method that is high-order accurate for smooth solutions and satisfies 



a discrete local conservation law like (1.3) and a discrete local positivity condition 
Ik — ^ "^ticre K is restricted to a union of mesh cells. 



Fig. 1.1: Assumptions and requirements of positivity limiting (Section^^. 



Accuracy is not the only goal of a numerical method, however. One also seeks 
physicality of numerical solutions. If a physical condition is maintained by the so- 
lution, then it is desirable for a numerical method to mimic this by maintaining a 
discrete version of this condition. Such methods are referred to as mimetic or con- 
forming. 

For example, in a problem where material is conserved, the solution u is restricted 
for all time to a manifold of solutions that have the same total amount of material. 
A numerical method that fails to satisfy a discrete material conservation law can over 
time drift from this manifold, resulting in unphysical behavior. Although one can 
correct for this by a global adjustment to the solution, such an approach still can 
yield locally incorrect physics. Specifically, solutions that fail to satisfy a discrete 
local conservation law result in simulated shocks that travel at incorrect speeds (see 
e.g. pages 237-239 of LeVeque [H]). 

Similarly, in a problem where the solution should remain positive, a numerical 
solution that fails to maintain a discrete positivity condition could drift from the set 
of positive solutions, resulting in an unphysical or even unstable solution. Again, 
one could globally damp the deviation from the global average. But to ensure local 
physicality and stability, positivity preservation should be enforced in a local manner. 

Finite volume methods are called conservative because they are designed to mimic 
the conservation property of conservation laws. Positivity-preserving methods are 
designed to mimic the positivity-preserving property. The challenge of positivity 
limiting is to design finite volume methods that are conservative, positivity-preserving, 
high-order accurate, and numerically stable. 

Generalizing to hyperbolic systems, the fundamental assumptions and require- 
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ments for positivity limiting are summarized in Figure [lT] 

1.1. Historical overview. Consider a ID conservation law 



dtu{t,x) +d^f{u) = 

which maintains the condition m > 0. 

A finite volume method partitions the domain into intervals called mesh cells. On 
a mesh cell of width Ax centered at x = Xi we denote the numerical cell average at 
time r by Uf: 



- / U{r,x)dx. (1.4) 



U" : 
' A 

An Euler step for a method-of-lines finite volume method (e.g., DG or WENO) up- 
dates the cell average as follows: 

[7-1 = [7^ - |i [h{Ur^,,„ Ut,i^) - h{Ur_,i„ UU,^)] (1.5) 

where U^_^i2 and U^_-^j2 approximate solution values on the left and right of cell 
interface Xi_ij2, respectively, and the numerical flux h{U^ ,1}^^) is consistent with 
physical flux: h{u,u) = f{u). 

In our abstract setting, what we know about f is that it is defined so that phys- 
ical solutions maintain positivity. Therefore, to maintain positivity in the numerical 
scheme, the numerical flux is defined in terms of the solution to a physical problem. 
This insight lead to the first general positivity-preserving scheme: the Godunov 
scheme [5], which iterates the following: 

1. In each mesh cell replace the solution with its average value. 

2. Physically evolve the solution for a time step At. 

If At is sufficiently short, then the flux at each interface is given by the solution to 
a Riemann problem; this is guaranteed if AtX < Ax, where A is the sum of left-going 
and right-going signal speeds propagating from each interface in the physical solution; 



see Figure 2.1 In this case, the Godunov scheme can be implemented by equation 



( 1.5 ) if the numerical flux h{U , J7^) is chosen to be the interface flux of the Riemann 



problem with initial states {U ,U^) and the solution is assumed to be constant in 



each cell. A simpler choice of h which maintains positivity in update (1.5) is the HLL 
numerical flux function [B]. The HLL flux is defined to account for the transfer of 
material implied by the solution to the Riemann problem modified by averaging in 
an interval that contains the signals emanating from the interface. To show that HLL 
preserves positivity, Perthame and Shu [TT| used a modified Godunov scheme: before 
averaging in each mesh cell, the physically evolved solution is first averaged in a set of 
nonoverlapping intervals each of which contains the signals emanating from one of the 
interfaces. The Godunov scheme is unfortunately only first-order accurate in space. 
To achieve high-order accuracy in space, the numerical fiux function h used in 



update (1.5 1 needs to accurately approximate the physical flux /. This will be the 
case if h is consistent with the physical flux and if U~ and are both high-order 
accurate approximations to the exact solution at the interface. Therefore, we assume 
in each cell a representation of the solution of the form 

U:{x):=U- + dU-, (1.6) 

where is the cell average and dJ7" is a high-order correction that is a polynomial 
of degree at most k. The DG method evolves such a representation, and WENO 
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Conservative local accuracy-preserving positivity limiting by 
linear damping of the deviation from the cell average 

U 
U 

u 

The positivity limiting framework introduced by Liu and Osher achieves high- 
order accuracy while maintaining positivity of the cell average U in each mesh cell. 
Prior to each time step the deviation from the cell average is linearly damped until 
a cell positivity condition is satisfied; a time step is then calculated that preserves 
positivity of the cell average. Zhang and Shu |16| reduce the cell positivity condition 
from positivity at every point in the mesh cell (i.e. for an infinite collection of point 
evaluation functionals) to positivity at a finite collection of positivity points. The 
limited solution is U := U + 9{U - U), where 6* e [0, 1] is just small enough so that U 
is positive at each (black) interior point and at each (gray) boundary node. 

Fig. 1.2: Positivity limiting illustrated for ID mesh cells with cubic polynomials. 




reconstructs such a representation with each time step. The flux function h remains 
unchanged (i.e. is defined using the Riemann problem with states and .) 

High-order accurate positivity limiters were developed in the seminal works of 
Zhang [15] and Zhang and Shu [16] . An in-depth review of the history of this method 



can be found in the excellent review article [18]. As depicted in Figure 1.2 in each 
mesh cell they linearly damp the high-order corrections just enough to enforce posi- 
tivity at a set of strategically chosen positivity points. For concreteness, we describe 
their method in detail for the ID scalar case. Let i^T = [xi - /S.xl2,Xi + Aa:/2] denote 
a mesh cell, and let {sj }"^q denote the points of the Gauss-Lobatto quadrature rule 
that has just enough points to exactly integrate polynomials of degree at most k. 



As depicted in figure 1.2 to limit i7", hnearly damp di7" just enough to ensure 



that the solution is positive at each quadrature point. In formulas: 

f/„u„:=mint/r(s,), e:=uAn\l 1, U: {x) := + 6 AU^ . (1.7) 

We note that B 6 [0, 1], and in particular, if J/min ^ 0, then = \. Finally, the cell 
averages are updated by using the damped solution in (|1.5|: 



c^r' = f/r-s^ [Mc^i;i/2>t^ii/2)-/^(c/^-i/2,^ti/2)], a-s) 

where Uf -^^j^ are the edge values computed from the limited solution \\.7\ . We will see 
that this algorithm maintains accuracy and guarantees the same positivity-preserving 
time step that is guaranteed if positivity is enforced everywhere in the mesh cell. 

Zhang and Shu have extended their positivity limiters to shallow water and gas 
dynamics |17j . as well as to higher-dimensional problems on both Cartesian and un- 
structured grids [T71 [2D] . 

As leveraged in [16) . a framework for positivity limiting can also be used to im- 
plement a maximum-principle-satisfying scheme. Suppose that the exact (entropy) 
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Outflow capping 

Cf7"+^ For an Euler step, the updated cell con- 

tent CU""^^ depends linearly on the time 
step At. A stable time step that caps 
outflow at e.g. a = 70% can be di- 
rectly calculated, and in the scalar case 
is Aisafo = min{^, Atstabic}- 



Fig. 1.3: Outflow capping in the scalar case. 




solution u{t, x) remains bounded in the interval [C/minj t^nax]- By negation and trans- 
lation of state space, one can use scalar positivity limiting to enforce that the affine 
functionals [/" - C/min and C/max - are both positive. 

1.2. Outflow positivity limiting framework. The main contribution of this 
work is an analysis and generalization of Zhang and Shu's positivity limiter based on 
limiting the potential and actual outflow from each cell. This yields a multi-stage 
framework we call outflow positivity limiting. 

1. The first stage is point-wise positivity limiting, which consists of enforcing 
positivity at a set of positivity points consisting of boundary nodes and inte- 
rior points. Enforcing positivity at boundary nodes ensures that numerical 
fluxes are defined with positive values. For shallow water and gas dynamics, 
when positivity of the depth or density is enforced at a boundary node, it is 
important to calculate fluxes with remapped states in order to desingularize 
wave speeds. Enforcing positivity at interior points may improve stability. 

2. The second stage is boundary average limiting. Like Zhang and Shu, 
we maintain positivity of the cell average by after each time step linearly 
damping the deviation from the cell average just enough so that a cell posi- 
tivity condition is satisfied. The essential difference is that our cell positivity 
condition directly caps the boundary average rather than requiring positiv- 
ity at points in the cell interior. Specifically, the optimal version of our cell 
positivity condition requires the boundary average to be no greater than the 
maximum boundary average possible if the solution were everywhere positive 
and had the same cell average. Enforcing positivity at the optimal interior 
points of Zhang and Shu is a means of boundary average limiting. 

3. The third stage is outflow capping. Rather than restricting the time step 
to the guaranteed positivity-preserving time step (and assuming that the 
estimated upper bound on wave speeds is legitimate), we directly calculate a 
maximal stable time step that caps cell outflow at e.g. 70%. The updated cell 
average varies linearly with the length of an Euler time step, and therefore 
outflow capping can be performed with the same limiting procedure that 
Zhang and Shu use to enforce positivity at a positivity point. See Figure [T3j 

1.3. Benefits of the outflow positivity limiting framework. This frame- 
work offers simplicity, insight, flexibility, and computational efficiency. Our framework 
is implied by a set of natural requirements, e.g. to maintain positivity of the cell aver- 
age and to guarantee a stable, positivity-preserving time step. We thereby decouple 
requirements and arrive at a framework that is intrinsic to the requirements and ar- 
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guably optimal. Outflow capping ensures that knowledge of the fastest wave speed 
is not needed to compute the positivity-preserving time step. A minimum positivity- 
preserving time step is guaranteed by positivity of the retentional, which is defined 
in terms of the data that actually determines the evolution of the cell average (i.e. 
boundary node data and the initial cell average), and has a physical meaning (cell 
content minus maximum possible loss) . The proof that positivity is maintained is thus 
intuitively obvious and works for spatially varying flux functions and for arbitrary cell 
geometry and representation space. 

Boundary average limiting can be used to guarantee a positivity-preserving time 
step even if the data used at the boundary nodes is not interpolated by a polynomial. 
This is relevant in the gas-dynamics case, where it becomes necessary to modify states 
at boundary nodes before evaluating fluxes, e.g. to desingularize fluid velocity when 
enforcing positivity of the density or to enforce realizability or hyperbolicity in higher- 
moment gas-dynamic models. 

Capping the boundary average is enforced via positivity of a single linear func- 
tional, the retentional, and is therefore no more expensive than enforcing positivity 
at a single point. For stabilization purposes, one may additionally enforce positivity 
at interior points, but the choice of these points can be optimized for requirements of 
efficiency, simplicity, or stability; if positivity points are instead optimized for max- 
imizing the guaranteed positivity-preserving time step, then enforcing positivity at 
these points enforces the optimal cap on the boundary average (as illustrated in Fig- 
ure 1.4 ) While we show that such optimal positivity points always exist (see Figure 
5.1), precisely isolating them can be non-trivial (in general requiring the solution of 
linear programming problems on a convex domain) , and in the case of nodal DG it is 
inefficient to check positivity at these non-nodal points if they are relatively numerous. 
In contrast, outflow positivity limiting merely requires an estimate of the maximum 
boundary crowding M* realizable by a positive solution. Note that precisely isolating 
M* also requires solving linear programming problems on a convex domain. 

1.4. Parts. This work is conceived of as the first in a three-part series. This first 
part is analytical rather than computational and avoids making claims dependent on 
computational experiment; rather, we lay out a general framework that incorporates 
previous work and that will guide subsesquent computational investigation. The 
second part will give a detailed justification of the analytical results. The third part 
is planned to consider issues that require computational simulation. In particular, the 
stability consequences of enforcing positivity at Zhang and Shu's points, at all nodal 
points, or only at boundary nodes, will be a natural follow-up study, and results will 
likely be dependent on the details of the choice of oscillation-suppressing limiters. 
Since our framework generalizes that of Zhang and Shu, it has in a sense already been 
computationally tested by their work [l4l ITSl fT6 l [TT l HS l [T9 l l20] . 

1.5. Section summary. The outfiow positivity limiting framework is developed 
in subsequent sections of this paper. In Section [5] we describe how DG (or WENO) 
updates the cell average. In Section [3] we show in the context of scalar conservation 
laws that enforcing positivity of the retentional guarantees that the cell average re- 
mains positive for a time step inversely proportional to the boundary crowding cap 
M chosen in the definition of the retentional. In Section IH we show that the same 
statement holds for hyperbolic systems of conservation laws. In Section [s] we develop 
a theory of admissible (accuracy-preserving) and optimal cell positivity conditions 
that shows that, as long as the boundary crowding cap is admissible (i.e. no less than 
an optimal M*), then enforcing positivity of the retentional by linearly damping the 
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Positivity limiting caps the boundary average. 




M*U 



BU* 




Retentional positivity limiting caps the boundary average by the maximum boundary 
average attainable by a positive solution. The original solution U = U + dU has the 
same cell average U as the retentional-positivity-limited solution U-jz = U + OndU and 
the pointwise-positivity-limited solution Ux = U + OxdU have. The optimizer U* 
shown has the maximum boundary average BU* = M*U out of all positive quartics 
with cell average U (Theorem 5.7). The interior (black) positivity points X = {±\/l/5} 



are optimal and are thus zeros of U* (by Theorem 5.6). At each positivity point, Ux 
is positive. For U-ji, the boundary average equals (at most) BU* and the optimal 
retentional H^'* U := M*CU - BU, proportional to the sum (in general a weighted 



sum) of the values at the interior points, is (at least) zero. See Figure 5.11 



Fig. 1.4: Retentional positivity limiting for a one- dimensional mesh cell with a quar- 
tic polynomial basis (M* =6), compared with pointwise positivity limiting. 



deviation from the cell average preserves the order of accuracy and can be accom- 
plished by enforcing positivity at strategically chosen interior points. In Section [6] 
we apply this theory of cell positivity conditions to compute optimal or near-optimal 
boundary crowding caps and optimal interior points for typical cell geometries and 
polynomial orders. In Section [7] we discuss practical implementation issues regarding 
efficiency, robustness, and generalization, with specific reference to shallow water and 
gas dynamics. These issues include efficient and correct positivity checks, wave speed 
desingularization, stability benefits of interior positivity points, stable optimal time 
stepping, multistage and local time stepping, and application to isoparametric mesh 
cells. In Section [8] we summarize the key results of this work. Boxed figures present 
a largely self-contained summary of the work. 

2. High-order discontinuous Galerkin schemes. For simplicity, we dis- 
cuss positivity-preserving limiters with reference to the discontinuous Galerkin (DG) 
method. This does not entail loss of generality, since a WENO finite volume method 
can be thought of as a DG method that reconstructs the high-order component of 
the representation prior to each time step. Standard DG schemes discretize the weak 



form of hyperbolic conservation law (1.1) 

d_ 
dt 



/ uip= / f-yip- <£ n-f(z), (2.1) 

Jk Jk JdK 



where K is an open subset of M^, dK is the boundary of K, h is an outward pointing 
unit vector that at every point on dK is perpendicular to the boundary dK, and 
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The Riemann problem with states {U ,U^) and flux function / is defined to 
be the problem 

(PDE) dtu{t,x)+d^f{u) = 0, u{0,x) = 

Riemann problems (and thus their solutions) are by requirement self-similar, i.e. 
invariant under dilation of space-time: u{at,ax) = u{t,x) for any a > and t > 0. 
So we can write u(t,x) =■ u{x/t) for t > 0. In hyperbolic problems information 
propagates at finite speed. Define the signals s^{t) < < s^{t) emanating from the 
interface x = so that a := [s~,s^] is the smallest interval such that u{t,-) agrees 
with u(0,-) outside a. By self-similarity, the signal speeds are independent of t 
and define the right-going signal speed s^(C/ , C/^) and left-going signal speed 
s~{U~ , U"^). The interface flux fa is defined to be the flux at x = needed to account 
for the amount of material that u accumulates on each side of the interface; it is 
independent of t > and equals /(u(0)) unless u happens to be discontinuous at 0; 
then subtracting (PDE) from f^^ (PDE), applying the fundamental theorem of 

calculus, and solving reveals that 2/o = f{U^) + f{U^) + dt[ u{t,-) - u{t, ■)] if 
x~ < ts~ and x'^ > ts^ . 



Fig. 2.1: Definition of the Riemann problem used to define numerical flux 



{ U- ifa;<0, 
1 [/+ ifa;>0. 



(y9 6 Cg is a test function. 

Assume that Vl is the union of a finite set of non-overlapping mesh cells. We 
assume that the solution, when restricted to any mesh cell K e 0,^, belongs to a 
polynomial space V = V{K) and is otherwise unconstrained (e.g. by the requirement 
of continuity across mesh cell boundaries) . Often V is P|, , the set of polynomials in D 
variables of degree at most k. More generally, we define deg(V) > to be the largest 
k such that V contains P^. The discontinuous Galerkin method approximates the 



exact solution to (2.11 with a piecewise polynomial function: 



N(k) 

E (2-2) 

where, on each cell K, {'^}^Sq^ is a basis for V{K) and {ipj}^}^'^ is its co-basis. The 
basis and co-basis are mutually orthonormal in the inner product IxV^Vj - 
where \K\\s the measure of the mesh cell K. 

In nodal DG, the co-basis consists of a set of point-evaluation functionals that 
evaluate the solution at the nodes, and the solution is thus represented by nodal values. 
In modal DG, the basis is often chosen to be a sequence of orthonormal polynomials 
of increasing degree. As a consequence, the basis and co-basis are identical and the 
coefficients of higher-order basis functions decay: specifically, for a smooth function 
M, the L^-projection coefficients fj^uip^ decay like 0{lS.x^ ^^), where Ax is the mesh 
cell diameter and k' is largest such that P|, c spanjiys' : i < j}. 

The DG method solves a discretization of equation (|2.1| in each mesh cell K: 

d 
At 



f Uifi^ = ['^{{t,jc,U)-Vip^ - f'^ h{U-,U\h)ip^, (2.3) 

JK JK JdK 
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Framework of a method-of-lines finite volume Euler step 

Integrating the hyperbohc differential equation 

dtu{t,x.) + V-f(t,x, u) = 
over a mesh cell K gives the integral form 

— f u + £ fi'f = 0. 

dt Jk JdK 

Making the replacements u ^ U e V, B, and n-f -* h, where V is a finite- 

dimensional representation space, his a numerical flux, and 6 is a numerical boundary 
quadrature, gives the ordinary differential equation of a method-of-lines finite volume 
approximation dtCU = -Bh, where C := An Euler step for this ODE is 

C[/"+^ =CU- AtBh. 

Typically V when restricted to a mesh cell is a polynomial representation space con- 
taining all polynomials of degree at most k, h = ZxeQ Wx/i(x) is a numerical 
quadrature rule with points Q c dK and weights Wx > 0, and h{x) is a numerical 
interface flux for the Riemann problem with frozen flux function / := u i-^ n - f (i, x, u) 



and states {U',U^) (see Figure 2.1); here C/"(x) is the value of C/ at x approached 
from within K and C/^(x) is the value of C/ at x approached from outside K. DG 
and WENO are two ways of updating the deviation of U from the cell average. 



Fig. 2.2: The context of this work as outlined in SectioniB 



where n is the outward pointing unit normal to dK and where h{U~ , U'^ ,n) denotes 
a numerical flux function where represents the state at the boundary node ap- 
proached from inside the cell and represents the state at the same boundary node 



approached from outside the cell. See Figure 2.2 In equation (2.3) we have replaced 
exact integration by quadrature rules; in particular, denotes a quadrature rule 
exactly equal to for polynomials of degree at most 2k - I and denotes a 
quadrature rule exactly equal to for polynomials of degree at most 2k. If i^T is a 
polytope then the boundary quadrature rule can be written as the sum over all faces 
of a Gaussian quadrature rule on each face. 

An important property of the DG scheme is that it conserves the total amount of 



U from time step to time step; indeed, taking (p = 1 eV, equation (2.3) becomes the 
conservation law 



A f U=-(f^h. (2.4) 
dt Jk JdK 



The entire context of this work is concerned with maintaining positivity in a single 
arbitrary mesh cell. We therefore adopt a simplified notation that omits explicit 
reference to the mesh cell K. We define the cell content to be the cell integral 
C := U <-* f^U or its scale-invariant version, the cell average C U ^ fj^ U . We 
define the boundary sum to be the boundary integral quadrature B := U ^ <f^ U~ 
or its scale-invariant version, the boundary average quadrature B := U i-^ f^j^ . We 
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use V := \K\ to denote the "volume" of the mesh cell and A := \dK\ to denote the 
"area" of its boundary. With these conventions, the ODE (2.4) is expressed as 



dtCU=-Bh, i.e., 
dtCU = -Bh, 

where h:= yh is a, scale- invariant version of the numerical flux, and an explicit Euler 
step is expressed as 

CU''^^ = CU - AtBh, i.e., (2.5) 
CU''^^ =CU-AtBh. (2.6) 



This method-of-lines finite volume framework is summarized in Figure 2.2 



3. Framework for outflow positivity limiting. In this work positive means 
nonnegative unless qualified with the adjective strictly. A positive combination means 
a linear combination with positive coefficients, and a positive functional is defined to 
be a functional that is positive on positive functions. 

A simple algorithm that maintains positivity of the cell average is to repeat the 
following sequence: 

1. Assume that the cell average is positive. 

2. If necessary, linearly damp (rescale) the deviation from the cell average just 



enough so that a cell positivity condition is satisfied. (See Figure 1.2 ) 
3. Execute an Euler step for a stable time step that is just short enough so as 

to guarantee that positivity of the cell average is maintained (given that the 

cell positivity condition is satisfied). 
This algorithm was first introduced by Liu and Osher [10 and further developed by 
Zhang and Shu [T6] . 

Assuming that numerical fluxes are calculated by evaluating the flux at the bound- 
ary node states (rather than at remapped states) , the cell positivity condition should 
at least require the solution to be positive at the boundary nodes so that the Riemann 
problem used to define the numerical flux function is well-deflned. To ensure that pos- 
itivity limiting does not compromise order of accuracy, it will be enough to require 
that the cell positivity condition be satisfied if the initial data is already everywhere 
positive (as we will verify in Theorem |5.2[ ). 

Given the constraints of this framework, we ask: What cell positivity condition 
and time step will guarantee that U"^^ > 0? If we are not allowed to modify a 
positive solution, then the maximum time step for which we can hope to guarantee 
positivity is the maximum such time step if the initial data is positive. We therefore 
ask: What is the maximum time step for which it is guaranteed that [/"^^ > if 
U >0 in K? We can assume that the outfiow h is strictly positive, since otherwise 
the cell average will remain positive for any positive time step. In this case, setting 
CC/"+i > in equation (IZSl) or (IzB gives 



. 1 . 1 Bh , 

Ar^ > At:^,„ := — = ^ (3.1 
CU CU 

That is, the time Atzcm until strict positivity of the cell average is violated is the 
ratio of the total integral of U in the cell to the net rate of flux out of the boundary. 
Clearly, to guarantee positivity we need a constraint on the fiux out of the boundary. 
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A constraint that extends naturally to systems is to specify a cap on wave speeds. 
Abstractly we impose that 



h<XU-, (3.2) 

where we call A the speed cap; equivalently, h < XU^ , where A := is a scale- 
invariant version of the speed cap. This condition clearly holds for a scalar problem 
with a convex flux function (e.g. Burgers' equation) if A is the wave speed, and in 



Theorem 3.7 we show that for general systems A can be defined as a cap on the sum 
of incoming and outgoing signal speeds at each boundary node. 

Remark 3.1 (omitting bars). The form of the scale-invariant equations is 
identical to the form of the non-invariant equations, so we can choose whether to 
omit or retain bars when discussing general properties. 



3.1. Boundary crowding caps and the retentional. Equations (3.1) and 



(3.2) imply that positivity is maintained if 

Ar^ > xBu = XBU, (3.3) 

where we call BU := ^ or its scale-invariant version BU := ^ the boundary crowd- 
ing. Therefore, we can guarantee a minimum positivity-preserving time step by en- 
forcing bounds on A and BU. To maintain order of accuracy, these bounds need to 



be physically justified. Enforcing a cap on A is briefly considered in Section 7.2 The 
focus of this paper is on enforcing a justified cap on BU. 



Given that AI > BU, equation (3.3) says that positivity is maintained if 



At-' > At-^l := AM. (3.4) 

To determine a justified M, we use that physical solutions satisfy positivity: 

Definition 3.2 (M*). We define the positive solutions to be the set := 
{U e V ■ U > in K and CU > 0} of representations positive and somewhere nonzero 
in the mesh cell, and we define the optimal boundary crowding cap (or opti- 
mal interior weight) M* to be the maximum boundary crowding over all positive 
solutions: 

M^(V) := sup B{U). (3.5) 

We can enforce that M > BU in exactly the same way that positivity is enforced at 
positivity points if we define this condition in terms of positivity of a linear functional: 

Definition 3.3 (retentional). Enforcing a cap M on BU = is equivalent to 
enforcing positivity of the retentional^ TZ^^U := MCU - BU, or one of its reseated 
versions such as TZ^ := CU-WBU , where W '■= M . We call M the interior weight 
( or boundary crowding cap ) and we call W the boundary weight ( or positivity 



CFL cap — see Theorem 3.12) 



^ An abbreviated form of retention functional, in the admittedly hokey mathematical tradition 
of using "adjectivals" as substantives and thus nouns. 
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Convention 3.4 (W = M ). Throughout this paper, any symbol involving the 
letter "W" represents the reciprocal of the corresponding symbol obtained by replacing 
"W" with "M". 

Definition 3.5 (admissible weights). If M > M* , then we say that M and W 
are admissible weights, because in this case positivity of U implies positivity ofTZ'^^, 
which means that enforcing positivity of 7?,*^ respects accuracy. 

Formally, we have the following two theorems, which comprise the essence of the 
paper. 

3.2. First fundamental theorem of outflow rate limiting. 



The first essential theorem of outflow rate limiting is displayed in Figure 3.1 It 
asserts that at each boundary node the numerical outgoing flux rate h is bounded 
by the product of the interior value at the node and a speed cap A; it assumes 
that boundary node states are positive and that h preserves positivity for a one-cell 
problem with states (0, U~ ,U^) and cell width Aa; = A for any At < 1. 

Remark 3.8. Multiplying by ceh area over cell volume, h{U',U^) < AC/". 

Remark 3.9. Examples of such a positivity-preserving flux function h are the 
numerical flux defined by the exact Riemann solver or the Harten-Lax-van Leer (HLL) 
[6] or local Lax-Friedrichs (LLF) 13J approximate Riemann solvers used with numer- 
ical signal speeds that are truly upper bounds for physical signal speeds. See Figure 
|3.1[ LLF is the special case of HLL where the left-going and right-going numerical 
signal speeds are equal: S~ + S'^ = 0. 

Remark 3.10. If outflow capping is used, then it is not necessary to compute 
the speed cap A; it is sufficient to know that such a finite cap exists. 

Remark 3.11. The proof assumes that Riemann problems are well-defined for 
vacuum states and involve finite wave speeds. Riemann problems involving vacuum 
states were considered for gas dynamics in |9]. For the proof it is enough that there 
exist a sequence of strictly positive states approaching vacuum such that the 
limsup of the signal speeds of the Riemann problems with states {Zm, U^) is finite. In 
systems such as shallow water and gas dynamics for which any state can be connected 
to the vacuum state without shocks, one can choose states Zm so that s^ {Zm,U~) = 
|AUax(C/") for all Z„; that is, s+(0,[/-) = |A|„,ax(C/-). 

3.3. Second fundamental theorem of outflow rate limiting. 

The second essential theorem is displayed in Figure |3.2| and states that if the 
retentional 7?.^ = C - WB is positive then an Euler step maintains positivity of the 
cell average for any time step whose nondimensionalized version XAt is no greater 
than W, where A is a uniform upper bound on the speed cap over all nodes. 

The first fundamental theorem can be invoked to satisfy the first hypothesis {h < 
AC/") of the second fundamental theorem if at each boundary node the numerical flux 



function is defined using a ID Riemann problem with frozen flux (see Figure 2.2 ). The 
proof of Theorem |3.7| assumes that physical solutions to the two auxiliary Riemann 
problems are well-defined and maintain positivity. For simplicity, we can impose this 



assumption rather than assume that the full conservation law (1.2) of Figure [O] 



maintains positivity. But in fact, omitting careful justification, it is an implication of 
this work that, with appropriate regularity and convergence assumptions, we have: 
Proposition 3.13. The following are equivalent: 



1. The full system (1.2) maintains positivity. 

2. Any Riemann problem with positive states and frozen flux function maintains 
positivity. 
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I: h(U ,U*) < XU for a positivity-preserving numerical flux. 

U- -\ 1 1 -^^ZT" 1 u{0,x) 

u{At,x) 

MHLL(At, a:) 




— (u-)-ll 

Definition 3.6. Define the one- cell problem with states {Z,U^ ,U^) , cell width 
Ax, and flux function f by 



dtu + dxf{u) = 0, 



i{Q,x) 



if X < 0, 

if < a; < Ax, 

if Aa; < x: 



(3.6) 



it is implied that Z,U ,U^ > 0. Define the speed cap to be the sum of the signal 
speeds entering the ceU: A( Z, , {7+ ) := s+ ( Z, ) - ( L/" , C/+ ) . 

Theorem 3.7 (outflow rate is bounded by speed cap times interior value). Consider 
the one-cell problem with states {Z,U^ ,U^), cell width Ax, and positivity-preserving 
flux function f . Define the numerical speed cap X'-= s^{Z,U^)-S^{U^,U^), where 
{U' ,U^) < {U^ ,U^) < 0. Let ho{A,B) be the interface flux of the Riemann 
problem with states {A, B) and let h(A,B) be a consistent numerical flux function 
that preserves positivity for the Euler update 



{u-y 



At 

Ax 



[h{U-,U^)-ho{Z, [/-)]> 



(3.7) 



if XAt < Ax, i.e., if incoming signals do not cross. Assume that Z, ZX{Z, Z,U ), 
and f{Z) all equal or approach 0. Then h{U^ , U^) < AC/". 



Proof Choose At = Aa;/A. Then (|3J]) says that h{U-,U^) < XU' + ho{Z,U-). But 
since material cannot flow out of a vacuum, ho{Z,U') < 0; to verify this rigorously, 
make the replacements U~' ^ Z, ->■ U~ , h ->■ Hq, and X ^ X{Z, Z,U~), and use that 
hoiZ, Z) = f{Z), which approaches 0. Therefore, h{U- , < XU' . □ 



To see how to define A and h so that the Euler update ( |3.7[ ) maintains positivity, 
for < t < At, let 3^(1) < < SQ{t) be the signals emanating from a; = and let 
s"(t) < Aa; < s^{t) be the signals emanating from x = Ax (see Figure 2.1). For the 
HLL flux, the numerical signal speeds and are required to satisfy < s^ < 
< < . Assume that A is an upper bound on the sum of the signal speeds Sq 
and -S^ . Since Aa; = AAt, the signals Sq and S^{t) := Aa; + tS^ do not cross (i.e., 
s5(At) < 5" (At) < s"(At)). An exact Riemann solver uses the flux of the exact 
solution at the Ax interface. The HLL solution UHLL(i,a^) equals u{t,x) except in 
the interval S{t) ■■= [S^ (t), (t)], where it equals the average value of u in S{t), 
U*{U-,U^) ■■= f(un-f(up+S*u+-s-u- ^ rj,-^^ needed at a; = Aa; to account for the 
amount of material that uhll accumulates on each side of the interface defines the 
HLL numerical flux h{U- ,U^) := \{f{U^) + f{U-) + {S^ + S-)U* - S^U^ - S-U-) 
(an instance of the formula at the bottom of Figure [2^. 



Fig. 3.1: A signal speed cap and boundary node value bound the outgoing flux rate. 
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II: A speed cap and positive retentional guarantee a 
positivity-preserving time step. 

Theorem 3.12. Suppose that: 

h < XU^ (numerical flux is bounded via the speed cap X), 
A At < W (positivity CFL number is bounded by W ), and 
WBU < CU (the boundary crowding is capped by M = W^^ ). 

Then the updated cell average CU - AtBh is positive. 

Proof. AtBh < XAtBU < WBU <CU. a 

Remarks. The proof says that the loss is at most the cell content. The retentional 
TZ^U ■■= CU - WBU represents the content retained if the maximum possible loss 
occurs and is posit ive precisely when the boundary crowding BU /CU is capped by 



M. Theorem 3.12 is sharp: a smooth example with spatially varying flux function 
satisfying f = hXu on dK shows that loss can equal cell content. Enforcing positivity 
of TZ^'^ maintains accuracy for scalar positivity limiting if M > M* . 



Fig. 3.2: A speed cap and boundary average cap bound the outgoing flux. 



3. The one-cell Euler update ( |3.7[ ) of Figure 3.1 with an HLL numerical flux 
function maintains positivity for sufficiently large X. 
Justification. To see that [l] implies [2] assume without loss of generality that the 
flux function is frozen at to = and xg = and in the direction n. The self-similar 
solution u{t, x) = u{x/t) to the auxiliary Riemann problem with positive states {A, B) 
can be approximated arbitrarily well by the solution u(t, x) to a multidimensional 



problem (1.1) with initial data mo(x) equal to B if n-x > 0, else equal to A. In 
particular, u{£,) = limt^ou(t,t^n). We omit a careful justification of this limit, and 
simply note that / evaluated at fi • x approximates n • f arbitrarily well for sufficiently 
small t and x; here we rely on the fact that f is differentiable (hence continuous) and 
that u by definition depends continuously on f wherever u is well-defined. Since u 
satisfies positivity, so do u and u. 

That |2] implies |3] follows from the definition of HLL in terms of averagings of 



physical solutions (see Figure 3.1) 



That [3] implies [T] follows from the assumption that the positivity-preserving algo- 
rithm described in this work converges at least in the first-order case of a representa- 
tion space that is constant in each cell. We remark that condition [3] gives a practical 
means of verifying that a physical system maintains positivity. □ 

3.4. Direct proof for the ID case. See Figure [3731 

3.5. AfRne-invariant definitions of wave speed and retentional. We have 
used bars in Sections [2H3] to indicate the scale-invariant formulation. Note that all 
definitions and statements remain valid if bars are dropped. Alternatively, we can 
redefine barred quantities using affine-invariant definitions. 

Heretofore we have made no distinction between physical and canonical mesh cell 
coordinates. In general, each physical mesh cell is the image of a canonical mesh 
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Comparison with Zhang and Shu (ID case) 

A corollary of the theorems in Figures [O] and [3^ is the ID case: 



Corollary 3.14. The Euler update (1.8) maintains positivity of the cell average if 



the retentional TUJl^ := U^-WBU^^ is positive in each cell and if the time step satisfies 

^^1^ < W , where A/2 is an upper hound on signal speeds and BUi = '^^^^ is the 

boundary average; we assume that the frozen flux function /ii+i/2 used at any interface 
is positivity-preserving, in the sense that if it is used at all interfaces then it maintains 
positivity of cell averages for data constant in each cell if A-^ < 1 (i.e., Ai is short 
enough that signals from cell interfaces cannot cross in Godunov's method). 

To facilitate comparison with Theorem 2.1 of |17| . we offer a direct proof. 

Proof. Substituting J7" = TZUi + WBUi (a boundary- weighted quadrature rule), 

> HU + ^ ^"1/2 ~ [^»+l/2(^i+l/2' t^i*fl/2) ~ '*i+l/2(2'+, 

^ [+C^tl/2 - ^ [Vl/2((7tl/2,^-) - Vl/2(f/^-l/2,C/tl/2)]. ' 

which is a positive combination of positive quantities if the Euler steps in brackets 
satisfy — ^* < 1; here we take Z+ and to be the vacuum state and we use that 
material cannot flow out of a vacuum: hi+ii2{Z+, U^^^j^) < and /ii-i/2(t^i^iy'2' •^-) - 
Zhang and Shu instead cLSSunie a, spa-tia-lly invRriciiit flux function i^h^j^ 

1/2 = ^i-1/2) and 

choose Z+ = U^_^^^ and Z_ = U^^^^^- They write the retentional as a weighted sum 



over interior positivity points (see Theorem 5.6 1 and choose W = W* (see Theorem 



6.2) 



Fig. 3.3: Direct proof for the ID case. 



cell under a coordinate map. In the case of isoparametric mesh cells, the coordinate 
map is a diffeomorphism, and it is important to work in canonical coordinates. If the 
coordinate map is affine, however, then we can avoid making a distinction between 
physical and canonical coordinates if we adopt an affine-invariant formulation. 

In canonical coordinates, the mesh cell is almost always a cube or simplex and 
the representation space almost always consists of polynomials. This remains true 
for isoparametric mesh cells. With an affine-invariant formulation, our results are 
generally applicable, independent of whether the canonical simplex is defined to be 
regular or the corner of a box. 

In Section [Hj we tabulate values or estimates of M* for regular canonical poly- 
topes. Thus, when using our tabulated weights to apply the outflow positivity limiting 
framework, one should take wave speeds to be in terms of their values in the coor- 
dinates of a regular canonical mesh cell. A natural way to do this is to use affine 
invariants. 



An affine-invariant formulation of the Euler step (2.5) that agrees with the scale- 
invariant formulation (2.6) in the case of a regular canonical mesh cell is CU"^^^ = 
CU - AtB^h^, where B^is an affine-invariant version of the boundary average that 
computes the arithmetic average over all faces of the average on each face; here /if := 
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' — is ail affine-invariant version of the numerical flux, where d^e is the area of 
the face of boundary node e and \T\ is the number of faces of the mesh ceh. Similarly, 
since wave speeds scale like fluxes, if Ae is a wave speed at boundary node e in the 
canonical coordinates of mesh cell K, then the quantity Ag := ^-^^y^Xe is the afflne- 
invariant version that agrees with the scale- invariant version Ae = -^Ag in the case of 
a canonical regular polytope. In the case of a regular canonical polytope the afflne- 
invariant quantities h^, Ag, and agree with their scale invariant equivalents hf,, 
Ae, and B, and the theory of Sections [3[|5| goes through. 

If one is computing with a non-regular canonical polytope K (e.g. the standard 
orthogonal simplex defined as the corner of a box) then an affine-invariant definition 
should be used in the definition of the retentional: 



TZ^^{U):=CU -WB^U. (3.8) 

If outflow capping is used, then computing wave speeds is not necessary, and use of 
the boundary face average B^ in the definition (3.8 1 of the retentional is the only 



modification needed to implement affine-invariant positivity-limiting. 

4. Systems versus scalar case. In the systems case of Figure the set of 
positive states 7-" is a convex set V c . A corollary of Theorem 3.4 in [12] is that 
any open or closed convex set is an intersection of half-planes. But any half-plane is 
the set on which some affine functional A is positive. Therefore, we assume that there 
exists a set V* of affine functionals such that a state u e is positive if for all A eV* 
A{u) > 0. We call A a state positivity functional. Any such affine functional A 
decomposes uniquely as A =: s + A, where s = A{0) is a scalar shift and A:= A - s is 
its linear part. Applying A reduces the systems case to the scalar case: 



5tA(u) + V-(Af) = 0. 

After applying A to states and A to fluxes, all statements and reasoning of Theorems 
|3.7| and |3.12| remain valid. For example, the inequality h < AC/" becomes the inequality 
(Ah) < Ayl(U"). 

Let 7?. : V ^ M be a retentional. Applied component-wise, the retentional defines 
a state- valued linear map TZ-V^ , which we identify with TZ. Applied point-wise, 
a state positivity functional A = s + A defines a map A = s + A : ^- V that we 
identify with A. Observe that TZA = ATZ. It is desirable that the retentional commute 
with all state positivity functionals: TZA = ATZ; then enforcing positivity of TZ is the 
same as enforcing TZA > for all A e T'* . This will hold if TZs = s, which holds if 
s is always (that is, if P is a convex cone) or if the retentional is rescaled so that 
TZl = 1. Therefore, for any retentional TZ we define its unital retentional to be 
TZ:= Then TZl = 1. To be concrete: TZ^^ = ^l^I^ ', we remark that TZ is an affine 
combination of state values and therefore is invariant under translation of state space. 

For typical systems such as shallow water and gas dynamics, the set of positive 
states is a convex cone (i.e. invariant under rescaling), so s = and A = A and there is 
no need to rescale the retentional to its unital version. We remark that by adding the 
trivial equation 9tUcxtra = to the system, where Ucxtra is a scalar taken to be 1 for 
physical solutions, defining the set of positive states to be {{ru,r) ■■ueT',r > 0}, and 
extending each positivity functional A= s + A to the linear functional A := (u, Mextra) 
■suoxtra + Au, the set of positive states can be assumed without loss of generality to be 
a convex cone. 
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State positivity limiting for gas dynamics 



energy (per mass) c>, 




The damping coefBcicnt O-p is defined so 
that Mp := M + 9-pAu 6 dV. There exists 
a state positivity functional Ap which is 
zero at u^. Here we take Ap = f^^, the 
energy in the reference frame of Vp. 




The damping coefficient O-p can be calcu- 
lated by finding a root of the quadratic 
polynomial pp |17j . By concavity of the 
pressure, positivity is also assured by us- 
ing the secant rule (i.e., by pretending 
that pressure is hnear); see [Tl]. Note 
that, before enforcing p > 0, enforcing 
p > may be needed [T7] . 



Fig. 4- 1- Depiction of positivity limiting of a state value for Euler gas dynamics. 
Positivity limiting requires determining where a linear path in state space intersects 
the boundary of positive states. 



While we can assume the scalar case without loss of generality when it comes to 
the positivity-preserving results of Section [3j the accuracy result of Theorem |5.2| does 



not go through e.g. to the gas-dynamics case (see Section 4.2 ) 



4.1. State positivity indicators. 

Outflow positivity requires the solution to the following problems; 

1. To cap outflow, determine the largest value of 6* e [0,1] for which u{d) := 
u + 9du satisfies positivity, where u = CU and du = lS.t^g_yJ3h. 

2. To cap boundary crowding at M, determine the largest value of e [0, 1] 
for which u{6) := u + 9du satisfies positivity, where u = TZU, du = TZdU, 
H := (^^1)-!^^, and := MC - B. 

3. To enforce positivity at Xg, determine the largest value of 6* € [0, 1] for which 
u{9) •■= u + 9du satisfies positivity, where u=U and du = d[/(xo). 

All these problems seek the intersection of a line segment with the boundary of a 
convex set of states designated positive and require computing a damping coefficient: 

Definition 4.1. Let V he a closed convex set of states deemed positive, and let 
V* he a collection of affine functionals such that V = {u^ : (V^ e V*) Au > O}. 
LetAeV*. Letu,ueM.^, where u^V. Define du:= u-u. Define 9 A{u,du) to he the 
largest 9 € [0, 1] such that A{9u + 9u) > 0, where we define 9 := 1-9. Define 9-p{u, u) 
to he the largest 9 e [0, 1] such that 9u + 9u e V. Ohserve that 9a{u,u) equals 1 if 
A(u) > and else equals ^i^"^^ = z^^- Observe that 9^^ (u, du) = min^^-p* 0yi(w, du). 

If T'* is a finite set, then it can be used directly to calculate 9-p, as illustrated 



in Figure 1.3 for the case of scalar outflow capping. In the example of Euler gas 
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dynamics, however, the state is positive if the energy in an arbitrary reference frame 
(a linear functional of conserved variables) is positive. But such energies comprise an 
infinite collection V* of linear functionals. Instead, one can use the density p (a linear 
functional), and pressure p (a concave nonlinear functional, see section 3.1 of |14] ) or 
pp (a quadratic functional of conserved variables, see [17 ) as a finite set of positivity 
indicators. See Figure ITT] 



4.2. Accuracy of positivity limiting in the systems case. For the case of 
scalar conservation laws, accuracy of pointwise positivity limiting via linear damping 
is established in Section [5] based on Theorem 15.21 

In the case of shallow water, positive states are characterized by positivity of a 
single linear functional (the depth). Therefore, the proof of accuracy for the scalar 
case can be invoked to conclude that positivity limiting does not compromise accuracy 
of the water depth. Assuming that fluid velocities are bounded, this in turn implies 
that positivity limiting does not compromise accuracy of the solution as a whole. 

But in the case of gas dynamics, one must also enforce positivity of the pressure. 
There exist physical solutions with arbitrarily large and rapid variation in density 
for which the pressure remains arbitrarily close to zero. Therefore, if positivity is 
enforced by damping the deviation of the cell average by a scalar value oi 9 e [0, 1], 
then a given amount of damping needed to enforce positivity of the pressure can entail 
damping of the density variation that is arbitrarily large in magnitude. 

On the other hand, if the pressure of the exact solution is strictly bounded away 
from zero, then positivity limiting respects accuracy for the simple reason that posi- 
tivity limiters are not triggered if the mesh is sufficiently fine, if the exact solution is 
smooth, and if M > 1. In this case, positivity limiting is really about practical robust- 
ness. The study of accuracy and stability for gas-dynamic problems (or sequences of 
problems) for which the pressure is not bounded away from zero is delicate and is left 
to future work. 

5. Admissible weights and cell positivity functionals. In this section we 
show that positivity of the retentional Ti'^^{U) can be enforced without loss of accuracy 
as long as the interior weight M is greater than or equal to an optimal value and 
that this holds precisely when positivity of TZ'^^{U) can be enforced by enforcing 



positivity at appropriately chosen interior points in the mesh cell K (see Figure 1.4). 

In Section |6] we use this theory (and interior points in particular) to compute (or 
isolate) the optimal value M* for important mesh cell geometries and representation 
spaces. As elsewhere, in this section the following context is assumed: 

Context 5.1. K is a compact mesh cell (in canonical coordinates) and V is a 
finite- dimensional polynomial space which contains the set P|) of all polynomials in 
D variables of degree at most k> 0. We always assume that U eV. 



Recall from equation (3.5 1 that the optimal interior weight M* = M^(V) is defined 



as the supremum of the boundary crowding B{U) := over the set of all nonzero 
solutions U e V positive in the mesh cell K , where C{U) := U and B{U) := /g^, U . 
The retentional TZwiU) := C{U) - WB{U), or U^' {U) := MC{U) - B{U), is positive 
for all such positive U precisely when M := W^^ is admissible, i.e., when M > M* . 

The following results will be proved in detail in Part II |7|. We summarize these 
results in Figure |5.1| and illustrate them in Figure |1.4[ 

Theorem 5.2 (An admissible weight preserves accuracy). Let M > M* . Then 
linearly damping the deviation from the cell average just enough to enforce positivity 
of the retentionalTZ^^U) retains order-k local accuracy. 
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Proof (sketch). Since physical solutions are positive, damping the variation from 
the cell average just enough to enforce positivity everywhere in the mesh cell preserves 
accuracy; a formal proof uses that the maximum n,„ax •= dJ7 i-*- max(dU{K)) and the 
magnitude of the minimum nmin '•= dJ7 max(-dC/(if )) are (equivalent) asymmetric 
norms on the finite-dimensional linear space of polynomials dU e V whose cell average 
is zero (see [E]). If M is an admissible weight, then the retentional TZ^^U) is positive 
if U is positive everywhere in the mesh cell, so enforcing positivity of TZ^^ (U) also 
retains accuracy. □ 

Theorem 5.3 (Enforcing positivity of the points of a nodal scheme is sufficient 
to enforce a positive retentional for some M < oo). Let X be a set of points capable 
of representing the solution and for which the cell average can be represented as a 
strictly positive combination of point values. Then < oo, where := supS(V^) 
and := {?/ 6 V : C;/ = 1 and (Vx e X) C/(x) > 0}. 

Proof (sketch). If U is positive on X then CU > 0, so B{U) is continuous on 
Vx- Using that Jimax and rtmin are equivalent asymmetric norms, one can show that 

■= {U eV^ - CU = 1} is bounded and hence compact. Therefore, B{Vx) has a 
finite maximum. □ 

Theorem 5.4 (A boundary-weighted quadrature rule gives an upper bound for 
the optimal weight) . Suppose that the cell content is given by a boundary-weighted 
quadrature rule, C{U) = Ti-iU) + WB{U), where the boundary weight W is strictly pos- 
itive and7l{U) represents a quadrature rule Zxejf "^xC^C^) for the retentionalTZwiU) 
with positive weights Wx that is exact for U eV. Then M* < M := . 

Remark 5.5. Enforcing positivity at the points X enforces positivity of TZw', we 
call X a set of retentional (or interior) points for the weight M. 

Proof. Indeed, 71 = C - WB is positive for all U positive in K, so for any positive 
f/, W-^ > = B{U), so W-'^ is an upper bound on M* . □ 

Theorem 5.6 (An optimal boundary-weighted quadrature rule exists). An op- 
timal quadrature rule C{U) = Y,x'^xU{x) + WB{U) exists such that M* = W^^ . 

Proof (sketch). If X is a subset of K, then 7?.*^ is positive on V^^- if M > ■= 
sn'pB{Vx) thus (if A" is a finite set) is representable as a positive combination 
of values at a subset A c X of points which comprise linearly independent point 
evaluation functionals; a compactness argument extends this statement from finite A 
to X = K. □ 

Theorem 5.7. An optimizer U* e exists such that M* = B{U*). 

Proof (sketch). By definition, M* is defined to be the supremum of the boundary 
crowding B{U) := ^^jjj over all nonzero U positive in the mesh cell K. Thus, any U 

positive on K gives a lower bound B{U) on ill*. The question is whether an arg max 
of B exists. Since B is continuous and is compact (see the proof of Theorem 



5.2), an arg max U* of B exists. □ 
Thus, it is possible to determine the optimal weight M* and a set of optimal 
retentional points simply by guessing an optimizer and an optimal boundary- weighted 
quadrature rule and using them to confirm one another. We now identify properties 
that reduce the set of guesses that we have to consider. 

Theorem 5.8 (Invariance under isometrics can be required of retentional points 
and optimizer candidates). 

Proof (sketch). A typical canonical mesh cell has a group of isometrics; by 
averaging over an orbit, for any solution U there exists a U invariant under the 
action of the group of isometrics of the mesh cell with the same values of C, B, and 
B. Similar averaging shows that every boundary-weighted quadrature rule C{U) = 
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Optimal accuracy-respecting cell positivity functionals 

Assume Context |5.1| X is a compact mesh ceU, V is a finite-dimensional polynomial 
space which contains the set of all polynomials in D variables of degree at most 

k>0,CU := 4- U, BU := fg^ C/, BU := CU/BU, V^:={UeV: inf U{K) > 0, CC/ > O} , 

M* := supB(V^O: and TZ'^' := MC - B. Then the following results hold: 

Corollary 5.11 (existence of an invariant optimizer and of invariant re- 
tentional points which isolate the optimal weight). Linearly damping the 
deviation from the cell average just enough to enforce positivity of the retentional TZ'^^ 
retains order-k local accuracy as long as M is admissible, i.e., M > M* (Theorem 
\5.^ . Any solution U eV positive in K and any correct boundary-weighted quadrature 

rule C{U) = Exex Wxf/(x) + WB{U) give a bracket B{U) <M*< W'^ for the optimal 
weight (Theorem 5.4-). Furthermore, an optimizer U* exists (Theorem 5.1) and an 



optimal boundary-weighted quadrature rule exists (Theorem 5.6) such that the bracket 
is an equality: B{U*) = M* = . Enforcing positivity at t he s et X of retentional 
points enforces positivity of the retent ional TZ^ (see Remark 5.5). The zero-set ofU* 



is nonempty if M* + 1 (Theorem 5.10) and must contain X if X is optimal (by Theo- 
rem 5.6). All these statements continue to hold if the optimizer U* or the quadrature 



rule (i.e. the retentional points X and the quadrature weights w^) is required to be 



invariant under the action of the isometrics of the mesh cell (Theorem 5.8). 



Fig. 5.1: Summary of the results of Sectiom^ 



T.xex Wxf/(x) + WB{U) has an invariant version C{U) = Zxex WxC^(x) + WB{U). □ 
Lemma 5.9 (The boundary crowding of a positive combination of two functions 
with positive cell average is a convex combination of their boundary crowding values) . 
Let U,Ve V^. Let s > and t > 0. Then B{sU + tV) = aB{U) + (1 - a)B{V) for some 
ae (0,1). 

Proof B{sU + tV) = f^fff = aB(U) + (1 - a)B{V), where a = ^cU%^- ° 
Theorem 5.10 (A non-constant optimizer has a zero). An optimizer U* has a 
zero in K unless U = 1 is an optimizer (in which case M* = 1). 

Proof We have that M* = BU* . Let U^in '■= miTiU*(K). Assume that Unun > 0. 



Let U := U* - Un,in- Since «([/„„„) = 1, by Lemma |5J) B{U*) = B{U + U^in) = 

aE{U) + (1 - a), for some < a < 1. If B{U) < 1 then B{U*) < 1 = B{Umin), 

contradicting the maximality of B{U*). If B{U) > 1 then B{U*) < B{U), again 

contradicting the maximality of B{U*). Therefore, B{U*) = 1. But B{1) = 1, so 1 is 
an optimizer. 

□ 

6. Calculated weights (functional analysis calculations). The previous 
section developed a general theory of cell positivity functionals. In this section, we 
use the results of Section [5] summarized in Figure |5.1|to determine retentional points 
and lower and upper bounds for the optimal weight Af for the two standard mesh 
cell geometries (box and simplex) and for polynomial representation spaces of varying 
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order. The most important results of this section are summarized in Figure [O] 

Remark 6.1 (retaining bars when computing weights). In accordance with Re- 
mark |3.1[ we dropped bars in the previous section without affecting the vaUdity of 
the general theorems. In this section, however, we obtain concrete numerical values 
for the interior weight M; since we want these results to be independent of the scale 
of the canonical mesh cell, we retain bars in this section. 

6.1. Definitions and summary of results. Recall that M*j^{V) denotes the 
optimal interior weight for mesh cell K with representation space V. We tabulate 
values or upper and lower bounds on the optimal interior weights for standard canon- 
ical mesh cells: the unit interval [0,1], the unit square [0,1]^ the unit cube [0, l]^ 
an equilateral triangle A^, and a regular tetrahedron A-^. Each upper bound M 
constitutes an admissible weight and is calculating using a quadrature rule for the 
retentional TZ'^ whose quadrature points can be used as positivity points. 

Definition 6.2 (canonical mesh cells). In M^, define canonical mesh cells: 

[0,1]^= regular box, = sphere, 

= regular simplex, : "star-regular" polytope. 

We define a star-regular polytope to be any star-convex mesh cell that can be 
centered on the origin to have a constant value of n-x, where n is the outward unit 
normal and x is position on the cell boundary. We use "j^f^ as a generic designation 
for a star-regular polytope. Regular polytopes are star-regular. Since most often the 
representation space is P|,, we define 



For practical use, we summarize our calculation of interior weights in Figure 6.1 

6.2. Additional definitions for regular polytopes. To calculate quadrature 
rules that yield admissible boundary weights for higher-dimensional mesh cells, we 
try to reduce to a one-dimensional problem using the following symmetry framework. 

Definition 6.3 (notation and conventions). Let C = denote the cell average. 
With star-regular polytopes in mind, assume that the canonical mesh cell is centered 
on the origin. Define the radius r to be fi-x, which for a regular polytope is the 
distance from the origin to the closest point on the boundary of the mesh cell. We 
will assume unless stated otherwise (and without loss of generality) that the radius 
of -fC is 1. Let Br denote the average over the boundary of the mesh cell rK rescaled 
to have radius r. 

Proposition 6.4 (The cell average is a weighted average of boundary averages 
over rescaled cells) . Let K be a star-regular polytope in D -dimensional space. Then 



- Lr^'-'Brdr 

C = - = D 

r^-idr 



[\^-^Brdr. (6.1) 
Jo 



Furthermore, if the polynomial representation space V is a subset o/P^ then Br 
polynomial in r of degree at most k. 



Proof. Equation ( 6.1 1 equates weighted averages and is justified by observing that 
the thickness of the infinitesimal shell [r, r + dr] ■ i^T is dr and its area is proportional 
to r^^^. Since any polynomial is a sum of homogeneous polynomials, to justify the 
final statement it is enough to observe that the statement is correct for homogeneous 
polynomials; indeed, if U is homogeneous of degree k' then so is Br. □ 
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Tabulated results of Section |6] 

6.1.1. Boxes and even star-regular mesh cells. If the mesh cell K is even 
(i.e., there exists K e K such that K c K- K), then for a regular polygon the optimal 
interior weight for polynomials of degree at most k is bounded above by the optimal 
interior weight -^qd for a £)-dimensional sphere. Based on equations (6.4) and (6.7), 
we tabulate exact values or bounding intervals for the optimal interior weight for an 
interval (M^p ^j), a square {M^^ ^^2), and a cube (M^^ ^-|3). Note that for a tensor 

product polynomial space the optimal weight is M^p independent of the dimension 
D of the box (see Remark 6.14). Note also that for fc < 3, enforcing positivity at the 
cell center enforces positivity of the retentional for the optimal weight (see Theorems 



6.6 and 6.7). 
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6.1.2. Simplices and arbitrary star-regular mesh cells. If the mesh cell 
is not necessarily even, then for a regular polygon the optimal interior weight is still 
bounded above by M^^. Based on equations (6.5) and (6.8), as well as Remarli 



6.10 and Theorem |6.9[ we tabulate exact values or bounding intervals for optimal 
interior weights for the triangle (M^a) and tetrahedron (M^a). Note that to enforce 
positivity of the retentional for the optimal weight, for a quadratic representation 



space it is sufficient to enforce positivity at the cell center (see Theorem 6.6), and, in 
the case of a simplex, for a cubic representation space it is sufficient to also enforce 
positivity at the center of each face (see Theorem 6.4). When using these weights. 



take the boundary average as the arithmetic average of face averages (5 3.5) 
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Fig. 6.1: Essential results of Section\^ A hounding interval is given where the exact 
value is unknown. 



22 



6.3. Linear and quadratic representation spaces. 

Theorem 6.5 (linear representation space). Let K be a star-regular polytope 
and let V = P}^. Assume that the origin is the only point whose orbit under the 
isometries of K is a singleton. Then the optimal interior weight is M* = 1 and 
the positivity-preserving time step is guaranteed simply by enforcing positivity at the 
boundary nodes. 

Proof. If C/ is a linear function then B{U) = 1, as can be seen by averaging over 
the orbit of U under the group of isometries of the polytope. □ 

Theorem 6.6 (quadratic representation space). Let K be a star-regular polytope, 
and let V = Pfj. Assume that the origin is the only point whose orbit under the 
isometries of K is a singleton. Then the optimal time step is guaranteed simply by 
enforcing positivity at the boundary nodes and at the cell center. Enforcing positivity 
at the cell center is equivalent to enforcing positivity of the retentional TZ* = M*C - B, 
where the optimal interior weight is M* = ^^j^ ■ 



Proof. As seen from equation (6.1), for any homogeneous quadratic U, BU = 
If [/ is a constant-valued function, then BU = 1. So the optimizer cannot be 



constant. So by Theorem 5.10 its zero set must be nonempty. We can require the 
optimizer to be invariant under the isometries of the polytope. So it must be a positive 
rotationally invariant quadratic that has a zero at the origin, i.e. a homogeneous 
quadratic. The support of the optimal interior sum is therefore restricted to the 
origin. So the retentional is proportional to the value at the cell center. □ 

6.4. Cubic representation spaces. 

Theorem 6.7 (odd-degree representation spaces for boxes). For boxes, ^"^^^d = 
"^[0*1]"' and likewise for other even polytopes and for spheres. In particular, for 
a cubic representation space, M* = , and positivity of the interior sum can be 
enforced simply by enforcing positivity of the value at the origin. 

Proof. The isometries of an even polytope by definition include negation x i-^ -x. 
Therefore an optimizer U* 6 P^"^-'^ symmetric under the isometries of the polytope 
lacks odd-degree terms, so U* 6 P|^. □ 

Theorem 6.8 (cubic representation space for simplices: optimal retentional 
points). Assume thatV = T'%. For any simplex, optimal retentional points are located 
at the cell center and at the centers of the faces. 



Proof. By Figure 5.1 the set of optimal retentional points must be zeros of an 
optimizer and can be required to be invariant under the isometries of K. Since V is 
cubic, this means that optimal retentional points can exist only at the center of the 
cell and at the center of its faces. (Note that an optimizer candidate uniformly zero 
on the boundary could not have a boundary crowding exceeding zero.) □ 
Theorem 6.9 (cubic representation space: optimal weight for triangle and tetra- 
hedron). Assume that V = P|). For a triangle the optimal interior weight is = 
20/9, and the optimal retentional points are located at the cell center and at the centers 
of the edges. For a tetrahedron the optimal interior weight is M^g = 11/6, 

Remark 6.10. In [20], Zhang, Xia, and Shu construct a boundary- weighted 
quadrature rule for triangles. For representation spaces P|5' and P^^^, their rule has 



interior weight Af^^ -^j = M^^^J = (ggg equation (6.2)), showing for example 



that <3 = M^o -^] 



Proof (sketch). By Theorem 6.4 and Figure 5.1 we may demand an optimizer 



that is invariant under the isometries of K and that is zero at the center of the mesh 
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Key results: optimal retentional positivity points 






1 1^ = 
1 1 


= 1.6 
= 2,3 
1 












For quadratic and cubic polynomial rep- 
resentation spaces, the set of optimal in- 
terior points for a box is s imp ly the cell 
center (Theorems 6.6 and 6.7). 



For quadratic polynomials in a simplex, 
the set of optimal interior poii its is sim- 
ply the cell center (Theorem 6.6 con- 
trast with Figure 1 in [18 ). 




V = P2(x) • P2(y) V = ¥^{x) - V^iy) 

For tensor product polynomial spaces, 
the optimal interior points for a box are 
indicated in Remark 6.14| (see [T7] ) . 



For cubic polynomials in a simplex, the 
set of optimal "interior" points consists 
of the cell center and the centers of the 



faces (Theorem 6.4 1. 



Fig. 6.2: Optimal interior points for standard mesh cells and low-order polynomi- 
als. Black positivity points are retentional (i.e. "interior") points. Gray points (when 
shown) are boundary nodes and depend on one's choice of a correct boundary quadra- 
ture. The black points are determined by mesh cell geometry and the representation 
space and are thus independent of the choice of gray points. Enforcing positivity at the 
black points automatically enforces positivity of the optimal retentional TV := M*C-B 
and guarantees that positivity of the cell average is maintained for the same mini- 
mum time step that one would guarantee by enforcing positivity everywhere in the 



mesh cell (assuming that the same wave speed caps are enforced, see Section 7.2). 
Positivity is also enforced at boundary nodes, but for a different purpose: to ensure 
that Riemann problems are solved with physical states. For a linear representation 
space, the set of retentional points is empty, and it is enough to enforce positivity at 
the boundary nodes. 



cell and at the centers of its faces. Both for a triangle and for a tetrahedron, up to 
multiplication by a nonzero scalar there exists a unique such cubic polynomial. This 
polynomial is definite (i.e. positive after negating if necessary), and evaluating B for 
this polynomial yields the values ^ for a triangle and ^ for a tetrahedron. □ 



Results for quadratic and cubic polynomials are summarized in Figure 
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6.2 



6.5. High-order representation spaces. 

Theorem 6.11 (weights for ID interval). 



^^[0,1] = ^^[0,1] = ^ •= 2 ■ ^ ' 

Note that the optimal boundary-weighted quadrature is the Gauss-Lobatto quadrature 
with n interior points and end weight W . 

Proof. For the canonical mesh cell K = [0, 1], the Gauss-Lobatto quadrature 

rQ " 

2/ U = W(U(Q) + U(l)) + ywiU(x,) (6.3) 

with n interior points is exact for polynomials of degree at most 2n + l. The function 
U*{x) = n"=i(a; - is zero at all interior points and is of degree 2n. So for the 
representation spaces Pj^ and P^'^^, the quadrature rule is exact, U* is a positive 
function in the representation space, and U* is zero on the interior points. Therefore 
the hypotheses of Theorem |5 . 1 1 1 are satisfied and the conclusion follows. □ 

Remark 6.12. For ID mesh cells, Zhang and Shu enforce positivity at the Gauss- 
Lobatto quadrature points, which of course implies positivity of the retentional |17| . 

The previous theorem is a special case of the following theorem. 

Theorem 6.13 (optimal interior weight for spherical mesh cells). 

Furthermore, M^^ ^ M^o, with equality precisely when D = 1 or k < A. 

Remark 6.14. For boxes one can also construct a boundary- weighted quadrature 
rule by taking the tensor product of a Gauss-Lobatto quadrature rule with a quadra- 
ture rule used to integrate over a face and averaging over all D such quadrature rules, 
as is done in for the case D = 2. The resulting interior weight is M* = M^g 
independent of D. This weight and set of quadrature points is optimal for a tensor 
product polynomial space; indeed, invoking Corollary |5.11[ a confirming optimizer is 
■■= Ylfii U^g*^^{xi), where C^[o'*i] ^ unit-interval optimizer. But for P^, in case 
D > 1 and A: > 1, the boundary weight of this theorem is an improvement over 

Proof (sketch). Recall from (6.11 that C = D r^^^Brdr, and observe that, for 



boxes and spheres, Br is an even polynomial in r of degree less than or equal to the 
degree k of the representation space. So we can write Br = B{r^), where degB = [k/2\. 
We seek a quadrature rule that maximizes the weight on the point r = 1. Making the 
substitution t = 1 - 2r^ shows that C is given by an integral with a Jacobi weight 
function. 

In case k = Am (or k = Am + 1) we use a Gauss-Radau quadrature rule for Jacobi 
weight (GRJ) with m interior points, which is exact for polynomials of degree at most 
2m (see JT^). In the case of a spherical mesh cell, an optimizer that confirms that 
the boundary weight is maximal is U* (r^) = YliliiT^ - ^iY, where the are the 
interior quadrature points corresponding to the interior quadrature points ti of the 
GRJ quadrature rule. 

In case k = Am + 2 (or k = Am + 3) , we use a Gauss-Lobatto quadrature rule for 
Jacobi weight (GLJ) with m interior points, which is exact for polynomials of degree 
at most 2m + 1 (see [1]). In the case of a spherical mesh cell, an optimizer that confirms 
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that the boundary weight is maximal is U* (r^) = YVili{r^ - rf^, where the are 
the interior quadrature points corresponding to the interior quadrature points ti of 
the GLJ quadrature rule. □ 
Theorem 6.15 (admissible interior weight for any star-regular mesh cell). 



Remark 6.16. Equality holds only for the case M° d = 1. Recall from Remark 



6.10 



that Af^" < A/^V^ < i!i±iKli±^^ which agrees with estimate (6.5) for even- 



dimensional representation spaces and for odd-dimensional representation spaces is 



slightly better than estimate (6.5 1 



Proof (sketch). We want a quadrature rule with such weight on the boundary 
average for C = D r^'^Brdr. Unlike for boxes and spheres, we cannot assume that 
Br is an even polynomial. For we use a GRJ quadrature (see with n points to 
a get a quadrature rule for C, and for P^^^ we use a GLJ quadrature (see [1]) with 
n points to get a quadrature rule for C; in each case, the resulting end weight is the 



reciprocal of the interior weight appearing in (6.5). □ 



6.6. Lower bounds for interior weights. In the previous sections we com- 
puted exact values and upper bounds for optimal interior weights. In the case of 
higher-order representation spaces where we have merely computed upper bounds, we 
would also like to have a lower bound that suggests how much the estimate might be 
improved. In particular, we will obtain bounding intervals that prove that 

M^Qn=n{k^), Mfg^jo =f)(fc'), and Min = n{k^). (6.6) 



Recall that for ID intervals the optimal weight is given by equation (6.2). For 
boxes and for simplices, we now bootstrap from this result to obtain lower bounds for 
the interior weight. 

6.6.1. Boxes. First consider the case of the box [0,1]'^. It has 2D faces, each 
of which is a box of dimension D -1. Each face has 2{D - 1) sub-faces, each of which 
is a box of dimension D - 2. Suppose that U is defined on the x = 1 face of the box 
and has interior weight ^^n-i when restricted to this face. Then U can also be 

thought of as a function on [0, 1]^ and has boundary weight given by the recurrence 
relation 

as will be shown in detail in Part II of this work [71; we can take M^^"^^ ■■= Af^^ 
By induction on D, the recurrence relation (6.7) implies that M^^ = 0{k^), so by 



Theorem 



6.13 



6.6.2. Simplices. Let X be a regular simplex of dimension D. Suppose that U 
is a polynomial of degree at most k defined on the face F oi K not containing vertex 
and suppose that U is positive when restricted to F . Then U extends uniquely to 
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Prescriptions for positivity limiting 

Check positivity efficiently. One can use interval arithmetic to inexpensively con- 
firm cell positivity in the vast majority of mesh cells. For nodal DG, choose positivity 
points to be nodal points to avoid additional computational expense. See §7.1[ 

Pad inequalities. In practice, error in machine arithmetic makes it problematic to 
enforce exact inequalities such as p > or p > 0. Instead, enforce p > tp ov p > Cp for 
some small ep > and > 0. See [2]. 

Use positivity points to improve stability. Enforcing positivity at additional 
interior points may improve stability (likely depending on the choice of oscillation- 
suppressing limiters). Since values used in volume and boundary quadratures must 
be calculated (or available) anyway, positivity of these values can be efficiently en- 
forced, and enforcing positivity at the points guarantees that state values used in the 
numerical method are always positive. For modal DG, positivity at any finite set of 
points can be efficiently checked, and it makes sense to include the optimal retentional 
points if they are known. 

Use wave speed desingularization for systems, especially for shallow wa- 
ter. When enforcing positivity in shallow water, fluxes need to be calculated with 
remapped states in order to desingularize wave speeds. See remarks in Algorithm [T] 
and Algorithm [2] and Section 7.2 For gas dynamics, wave speed desingularization is 



needed when enforcing positivity of the density (although more often it is positivity 
of the pressure that needs to be enforced). 

Estimate an optimal time step after enforcing the cell positivity condition. 

The factor by which the time step must be shortened determines the expense of pos- 



itivity limiting (5 7.31. Enforcing the cell positivity condition guarantees a minimum 



time step Aipos that preserves positivity of the cell average; see equation (3.4). One 



can directly calculate the maximum one-stage time step A^zoro that maintains posi- 



tivity of the cell average (S4.1) in addition to calculating the maximum stable time 
step A^stabic to determine an optimal safe time step. For multistage and local time 
stepping, one can use Aigtabic, Aizdo, and optionally Atpos to maintain and iteratively 
adjust an estimate of a safe time step that is both stable and positivity-preserving. 



Fig. 7.1: Key points of Section^^ 



a polynomial that is homogeneous of degree k when expanded about V . This yields 
a recurrence relation for a lower bound M^'^ < Af^^, 



M 



k- 



■D+k\r^+D 



D 



1 + D 



(6.8) 



as will be shown in detail in Part II of this work [7|- By ind uction on D, recurrence 
relation ([6^ imphes that M^^^ = 0{P). Thus, by Theorem 



6.15 



Af; 



n{k^). 



7. Practical application of positivity limiting. We briefly consider the most 
prominent issues that arise when enforcing positivity. Detailed treatment of practical 
issues is deferred to Part II of this work [7j . We summarize in Figure |7.1[ 
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7.1. Efficient implementation of positivity checks. For almost all prob- 
lems, in the large majority of cells one can confirm that the solution is positive and 
will remain so for any stable time step by using interval arithmetic. Thus, one can 
enforce positivity with little additional computational expense per time step. For ex- 
ample, in the case of gas dynamics, one can quickly confirm positivity of the pressure 
in the vast majority of cells by checking that Pmin •= Pmin^min - rnmax/2 is positive, 
where pmin and Synin are lower bounds on density and energy in the cell and m^ax is 
an upper bound on momentum. 

The question is how to obtain a tuple of intervals [Umim '^Lmax] ^'^^ components 
of the conserved variables. In modal DG one can obtain lower and upper bounds that 
hold globally in the mesh cell by using lower and upper bounds on the polynomial 
basis functions (pj. This works well, because for smooth solutions the coefficients of 
higher-order modes decay rapidly as the mesh is refined. Nodal DG represents the 
solution using values at nodes; therefore, one can easily obtain an interval [Uj^i^, u^^ax] 
that applies to all nodal values. Thus, even for nodal DG, checking positivity at all 
nodes is usually no more expensive than checking positivity at a point. By Theorem 
|5.3| positivity at the nodes of a nodal DG scheme implies positivity of a retentional 
with weight W < W* , because the cell average is a positively weighted sum of values 
at the nodes (including boundary nodes) and so is a boundary-weighted quadrature. 

7.2. Wave speed desingularization is needed when limiting density. Let 

p denote depth in the shallow water case or density in the gas dynamics case. After 
enforcing positivity of p at boundary nodes, wave speeds need to be de-singularized. 
This can be done by modifying states used to compute fiuxes when the wave speed 
exceeds an estimate of the maximum wave speed in the physical solution. 

In both shallow water and gas dynamics, desingularization of momentum (i.e. 
fluid velocity) is necessary. Let Ucap be a fluid speed cap that, for a sufficiently refined 
mesh, is guaranteed to exceed the fluid speed that arises in the exact solution. Then 
an accuracy-respecting desingularization map of the momentum (or fluid velocity) is 
the rescaling u <- u- i?i(ucap/|u|), where Ri{x) = 1 if a: > 1 and equals x (or the spline 
x{2 - a;)) if < a; < 1 (see Section 4.4 of [T]). 

For shallow water, enforcing positivity means enforcing positivity of p and there- 
fore essentially always entails wave speed desingularization. For gas dynamics, in 
addition to desingularizing the fluid speed u, it may also be necessary to desingular- 
ize the sound speed c = n/t^, e.g. by capping the pseudo-temperature 9 := p/p. 

More careful study of wave speed desingularization is left to future work. 

7.3. Cost of positivity limiting. For most problems, for the vast majority of 



cell updates, a quick check (^ 7.1 1 will confirm that positivity limiting is unnecessary. 
If local time stepping is used, then the total cost of positivity limiting will be marginal. 
But if global time stepping is used, then the expense of positivity limiting is measured 
by the factor by which the time step must be shortened to maintain positivity. This 
suggests (1) a careful consideration of how to obtain tight and reliable wave speed 
bounds for use in wave speed desingularization and (2) a comparison of stable and 
positivity-preserving time steps. 

For (2), in the case of one-dimensional mesh cells with fcth order polynomial space, 
values are known: an approximate value for the maximum stable time step is given by 
^^stabic ~ {k + 1/2) X/ Ax if an SSP-RK time-stepping method of order k is used (see 
Table 2.2 in [5]); in comparison, capping boundary crowding by M* = {n + l)(n + 2)/2 
gives AtpQg = AM* = {n + l){n + 2)X/Ax, where k equals 2n or 2n + 1. 
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7.4. Isoparametric mesh cells. For a typical isoparametric mesh cell, in canon- 
ical coordinates the mesh cell is a box or simplex and the representation space is a 
polynomial space (e.g. P|,). Each physical mesh cell is the image of such a canoni- 
cal mesh cell under a diffeomorphism 0, and the physical representation space is the 
push-forward under (j) of the canonical representation space. The outflow positiv- 
ity limiting framework is defined in canonical coordinates and can be applied without 
modification for isoparametric mesh cells: under reasonable regularity assumptions on 
(/), direct application of positivity limiters in canonical coordinates respects accuracy 
and guarantees a minimum positivity-preserving time step [7]. 

8. Conclusion. This work consists of two main results: a framework and algo- 
rithm for positivity limiting, and calculations of optimal weights and positivity points 
needed by this framework. 

8.1. Outflow positivity limiting framework. We have developed a frame- 
work for preserving positivity of each cell average based on limiting the rate and 
amount of material that can flow out of the cell. This is achieved by linear damping 
of high-order corrections, remapping boundary node states to limit wave speeds, and 
limiting the time step. In each cell, the high-order corrections are linearly damped 
just enough to enforce a physically justified cap M > 1 on the boundary crowding 



BA(U) 
CA(U) 

where C is the cell average and B is the arithmetic average over all faces of the average 



BA{U) := ^^^jj^ for all affine functionals A in a set of state positivity functionals V* 



over each face (see Section 3.8); this condition is satisfied precisely when the unital 
retentional 

_ M-CU-BU 
M-l 

satisfies positivity. One way to enforce positivity of the unital retentional is to enforce 
that the solution is positive at appropriately chosen quadrature points in the mesh 
cell, as Zhang and Shu have done [ini HH] ■ Directly enforcing positivity of the unital 
retentional is computationally inexpensive and makes it unnecessary to determine and 
enforce positivity at these points. For concreteness, we display positivity-preserving 
algorithms for scalar conservation laws and shallow water (Algorithm [T]) and for gas 
dynamics (Algorithm [2| . 

8.2. Optimal weights and positivity points. Just as point-wise positivity 
limiting prompts the search for optimal positivity points, retentional positivity limit- 
ing prompts the need for a reasonably close upper bound M on the optimal interior 
weig ht M^(V). In Section [H] we developed a general framework to estimate upper 
and lower bounds on AfJf(V). In Section [Hj we applied this framework to tabulate 
interior weights for boxes and simplices for polynomial representation spaces. 

In practice, the canonical mesh cell is almost always a cube or a simplex, and, for 
problems that are likely to entail positivity limiting, the representation space is likely 
to be a space of polynomials of at most cubic degree. Assuming these conditions, 
the set of optimal interior positivity points is very small, and the essential take- 
home result of this work is that to guarantee the same positivity-preserving time step 
as if positivity were enforced everywhere in the mesh cell, in addition to enforcing 
positivity (and, for shallow water, desingularizing wave speeds) at boundary nodes, 
it is sufficient to enforce positivity at the cell center, except that for a simplex with a 
cubic representation space one must also enforce positivity at the centers of the faces. 
See Figure |6.2[ 
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Fig. 8.1: Key definitions for cell positivity. 
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Algorithm 1 Positivity-preserving Euler time-step for scalar conservation laws and 

shallow water. 

For shallow water, U, h, and f represent the first component of U, h, and f, 
respectively. For the canonical mesh cell K and representation space V determine 
an admissible weight M such that Ti^'U ■■= MU - BC/ > VC/ = S^ U'ip^> (see 



Figure 6.11, where BU = 'ZjU^B{ipj) is the boundary average (see 5 3.5). For 



efficiency, precompute the values B{(pj). 

1. \/K rescale the deviation dU from the cell average U by a percentage 9y, damping 
just enough so that t7 > at a set of positivity points X that includes all quadrature 
points in and so that {U) = MU -BU>0: 

U:=U", U:=CU, dU:=U-U, C/„iin = min [/, 

xeX 

9x = minjl, -^-'^^ — |, 6i^n =min|l, ^^^^^^|, Oy ■= inm{ex,0^jj) , 



If X is sufficiently rich (see Figure 5.1) then 9x < ^7^11 and computing 6:j^m is 
unnecessary. 

2. Compute a safe time-step At: 

(At,ero)"' = max l fdKKU,U,n) \ ^ ^ max|(a,Ai,ero)"\ (Ai.table)"'} , 

all k\ JkIJ / ^ ' 

where Afgtabio is the maximum stable time-step and < < 1 is a safety factor. 

3. (for shallow water, not scalar case.) Modify the boundary node states Uq if 



necessary in order to desingularize wave speeds at boundary nodes (see 5 7.2) 
4. Vi^ update the solution: 

^_U"+V^ = ^U(/3^ + At^'^f(r,x,U)-V(p^ - At j^\(UQ,U^,n) 
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Algorithm 2 Positivity-preserving Euler time-step for gas dynamics (cf. §3.2 of [16]). 
Choose large enough values £p> and ep > that vanish with machine epsilon. 

1. VK evaluate the solution U" =: U = U + dU at a set of positivity points X 
which includes the set of boundary nodes Q (i.e. the quadrature points of /g^) to 
determine the nodal states ]J_^ d Uq . 

2. y K linearly damp the high-order corrections dU by a percentage 9p, damping 
just enough so that the density p is positive at the points X and so that the density 
retentional TZ^^{p) = Mp - Bp is positive: 

9x = inm\l,^ — —\, pmin = minp, 



p-p 



mm 



Jl_ if Bidp)<{M-l)-p-ep 

^^^M = |(M^ if Bidp)>{M-l)-p-ep, 

U ^ U + 6*^ dU, \lj^^]l+9p dV^ , 0p := mm{ex , O^j^j 



3. linearly damp the solution just enough so that the pressure is positive at 
the points X: 

0p := max {6* e [0, 1] : (Vx e X)p(U + 61 dU(x)) > ep} , 
U^U + 6'pdU, Uq ^U + 6lpdUQ. 

4. Modify the boundary node states Ug if necessary in order to desingularize wave 



speeds at boundary nodes (see S 7.2 1 



5. yK and for some admissible interior weight M linearly damp dU so that the 
retentional 7?,*^U := MU - f^]J_Q is positive: 

Oq := max 6 [0, 1] : p ((17 - 1)U - eiBdUg) > O} , 

U^U + ^gdU, Ug ^U + SgdUg. 

6. Compute a stable, positivity-preserving time-step At: 

At := a, max | At e [0, a^^ At.tabic] ^ P (^^^ U - At (Ug, U^) j > o| , 

where Atgtabio is the maximum stable time-step and < az < 1 is a safety factor. 

7. update the solution: 

= Uv3^ + At J^^ f (r,x, U) • Vifi^ -At£\ (Ug, U^, n) 



93-' . 



32 



